#include "field.h"

void field_divergence(Field       *f_dvg,
                      const Field *f_x,
                      const Field *f_y,
                      const double inv_dx,
                      const double inv_dy)
{
    CHECK_FIELD(f_dvg);
    CHECK_FIELD(f_x);
    CHECK_FIELD(f_y);
    CHECK_FIELDS(f_x,f_y);
    CHECK_FIELDS(f_x,f_dvg);
    
    {
        double **       d     = ARRAY_OF(f_dvg);
        const double ** Fx    = CONST_ARRAY_OF(f_x);
        const double ** Fy    = CONST_ARRAY_OF(f_y);
        const double    fac_x = 1.0 * inv_dx;
        const double    fac_y = 1.0 * inv_dy;
        const int       nx    = f_dvg->nx;
        const int       ny    = f_dvg->ny;
        const int       nx1   = nx-1;
        const int       ny1   = ny-1;
        int i,j;
        
        for(j=0;j<ny;++j)
        {
            const int jm = (j <= 0)   ? ny1 : j-1;
           // const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                const int im = (i <= 0)   ? nx1 : i-1;
             //   const int ip = (i >= nx1) ? 0   : i+1;
                d[j][i] =
                fac_x * ( Fx[j][i] - Fx[j][im]) +
                fac_y * ( Fy[j][i] - Fy[jm][i]);
            }
        }
        
    }
    
}
